#
# Copyright (c) 2009, 2010, Davide Cittaro
# All rights reserved.
# 
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#     * Redistributions of source code must retain the above copyright
#       notice, this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#     * Neither the name of the <organization> nor the
#       names of its contributors may be used to endorse or promote products
#       derived from this software without specific prior written permission.
# 
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#

import ConfigParser
import re, copy, sys
import numpy as np
from config import Enzymes, Modifications, Elements, Aminoacids, XLinkers
import helpers
from hashlib import md5
import scipy.stats as SS


class Interaction():
  """ Class aggregating multiple xlinks
  """
  def __init__(self):
    self._poschksum = ''
    self._GScore = 0.0
    self._xlList = np.array([])
    self._xPosInProt1 = 0
    self._xPosInProt2 = 0
    self._prot1 = None
    self._prot2 = None
    self._filtered = False
  def __cmp__(self, other):
#    return cmp(self.GScore(), other.GScore())
    return cmp(self.aggregateMatchScore(), other.aggregateMatchScore())
  def __eq__(self, other):
    if self._poschksum == other._poschksum: return True
    return False
  def addXLink(self, xl):
    """ add a xlink to this interaction
    """
    self._xlList = np.append(self._xlList, xl)
    self._prot1 = xl.lPeptide().motherProtein()
    self._prot2 = xl.rPeptide().motherProtein()
  def setXPos1(self, pos):
    """ set the position on protein 1
    """
    self._xPosInProt1 = pos
  def setXPos2(self, pos):
    """ set the position on protein 2
    """
    self._xPosInProt2 = pos
  def setPositionCheckSum(self, c):
    """ set the position checksum
    this will be the checksum for this interaction
    """
    self._poschksum = c
  def evalPositionCheckSum(self):
    """ calculate the position checksum as md5
    """
    lacc = xl.lPeptide().motherProtein().accession()
    racc = xl.rPeptide().motherProtein().accession()
    m = [lacc + ':' + str(self._xPosInProt1), racc + ':' + str(self._xPosInProt2)]
    m.sort()
    self._poschksum = md5(';'.join(m)).hexdigest()
    self.positionCheckSum()
  def positionCheckSum(self):
    """ return the position checksum
    """
    return self._poschksum
  def length(self):
    """ return the numner of matches for this interaction
    """
    l = 0
    for xl in self._xlList:
      l += len(xl.spectra())
    return l
  def linkedProteins(self):
    """ return the proteins involved
    """
    return [self._prot1, self._prot2]
  def xPositions(self):
    """ return the position for the interaction
    """
    return [self._xPosInProt1, self._xPosInProt2]
  def xlinks(self):
    """ return the xlinks. Note that each xlink
    may have multiple matches
    """
    return self._xlList
  def setGScore(self, p):
    """ set the score as raw probability value
    """
    self._GScore = p
#    self._score = -10 * np.log10(p)
  def GScore(self):
    """ return the score
    """
    return self._GScore
  def aggregateMatchScore(self):
    """ return the median of the intensity scores
    for all the matches in this interaction
    """
    xScores = {}
#    for xl in self.xlinks():
#      xScores[xl.checksum()] = np.mean(xl.matchScore())
#    return np.median(xScores.values())
    s = np.array([])
    for xl in self.xlinks():
      s = np.append(s, xl.matchScore())
    return np.median(s)
  def retainBestSpectra(self, options):
    """ for each match check that the p-value is
    in the "top list". options["retainBestSpectra"] models this:
    if the value for that option is n, then the method will check
    all the matches assigned to each spectrum for this match (i.e. 
    other xlinks) and get the n best scores
    """
    n = options["retainBestSpectra"]
    tempXList = np.array([])
    for idx, xl in enumerate(self.xlinks()):
      tempSpectra = np.array([])
      for sp in xl.spectra():
        sp.orderMatches()
        if xl.matchScore(sp) < options["mScoreFilter"] or xl.pValue(sp) >= options["pValueFilter"]:
          # don't need these
          sp.removeMatch(xl)
          continue
        scores = [x.pValue(sp) for x in sp.matches()]
        scores.sort(reverse=False)
        if n <= 0:
          threshold = scores[-1]
        elif len(scores) >= n:
          threshold = scores[n - 1]
        else:
          threshold = scores[-1]
        if xl.pValue(sp) <= threshold:
          tempSpectra = np.append(tempSpectra, sp)
        else:
          sp.removeMatch(xl)
#      if len(tempSpectra) == 0: continue
      # this is slow and not mem efficient
      # but it's the only way to make it work properly
      # until I get rid of it :-)
      txl = xl#.copy() 
      txl.setSpectra(tempSpectra)
      txl._MatchScores = {}
      txl.evalMatchScores()
      tempXList = np.append(tempXList, txl)
    self._xlList = tempXList
    if self.length() < options["interactionLengthFilter"]:
      self.filterOut()
  def isFiltered(self):
    return self._filtered
  def filterOut(self):
    self._filtered = True

class XLink():
  """Class storing basic xlink information
  """
  def __init__(self):
    self._lPeptide = Peptide()
    self._rPeptide = Peptide()
    self._mmass = 0.0
    self._xlinker = ''
    self._lProtein = Protein()
    self._rProtein = Protein()
    self._lMap = 0
    self._rMap = 0
    self._assignedSpectra = np.array([])
    self._xPosInPep1 = 0
    self._xPosInPep2 = 0
    self._pValues = {}
    self._MatchScores = {}
    self._checksum = ''
    self._poschksum = []
    self._weight = 1.0
  def copy(self):
    """ return a copy of a xlink having the same base
    properties"""
    cxl = XLink()
    cxl._lPeptide = self._lPeptide.copy()
    cxl._rPeptide = self._rPeptide.copy()
    cxl._mmass  = self._mmass 
    cxl._xlinker  = self._xlinker 
    cxl._lProtein  = self._lProtein 
    cxl._rProtein  = self._rProtein 
    cxl._lMap  = self._lMap 
    cxl._rMap  = self._rMap 
    cxl._assignedSpectra  = self._assignedSpectra.copy() 
    cxl._xPosInPep1  = self._xPosInPep1 
    cxl._xPosInPep2  = self._xPosInPep2 
    cxl._pValues  = self._pValues.copy()
    cxl._MatchScores  = self._MatchScores.copy()
    cxl._checksum = self._checksum
    cxl._poschksum = self._poschksum
    cxl._weight = self._weight
    return cxl
  def __cmp__(self, other):
    return cmp(self._lPeptide.sequence(), other._lPeptide.sequence())
  def __pcmp__(self, other, sp):
    return cmp(self._pValue[sp.checksum()], other._pValue[sp.checksum()])
  def setLPeptide(self, pep):
    """ set the first peptide (left)"""
    self._lPeptide = pep
  def setRPeptide(self, pep):
    """ set the second peptide (right)"""
    self._rPeptide = pep
  def setLProtein(self, prot):
    self._lProtein = prot
  def setRProtein(self, prot):
    self._rProtein = prot
  def lProtein(self):
    return self._lProtein
  def rProtein(self):
    return self._rProtein
  def setMass(self, mass):
    """ set the molecular mass
    for this xlink"""
    self._mmass = mass
  def calculateMass(self, options):
    """ calculate the molecular mass
    and set it"""
    self._mmass = self._lPeptide.molecularMass() + self._rPeptide.molecularMass() + XLinkers.getfloat(self._xlinker, 'mmass')
  def mass(self):
    """return the molecular mass"""
    return self._mmass
  def lPeptide(self):
    """return the first (left) peptide"""
    return self._lPeptide
  def rPeptide(self):
    """return the second (right) peptide"""
    return self._rPeptide
  def setXlinker(self, xlinker):
    """set the xlinkers for this xlink
    """
    self._xlinker = xlinker
  def xlinker(self):
    """return the xlinkers for this xlink
    """
    return self._xlinker
  def setLXPos(self, pos):
    """set the first xlink position
    relative to the first peptide sequence
    """
    self._xPosInPep1 = pos
  def setRXPos(self, pos):
    """set the second xlink position
    relative to the second peptide sequence
    """
    self._xPosInPep2 = pos
  def lXPos(self):
    """return the first xlink position
    relative to the peptide sequence
    """
    return self._xPosInPep1
  def rXPos(self):
    """return the second xlink position
    relative to the peptide sequence
    """
    return self._xPosInPep2
  def setSpectra(self, list):
    """set the spectra matching to this
    xlink molecular mass. It's a list
    """
    self._assignedSpectra = list
  def addSpectrum(self, sp):
    """append a spectrum to the assigned
    to this xlink
    """
    self._assignedSpectra = np.append(self._assignedSpectra, sp)
  def spectra(self):
    """return the assigned spectra
    """
    return self._assignedSpectra
  def numberOfMatchingSpectra(self):
    """return the number of spectra for this
    xlink
    """
    return len(self._assignedSpectra)
  def numberOfMatchingIons(self, sp = None):
    """Return the number of matching ions
    for all assinged spectra. If a spectrum is defined
    returns only the matching ions for that spectrum
    """
    if sp:
      return self.lPeptide().numberOfMatchingIons(sp) + self.rPeptide().numberOfMatchingIons(sp)
    return np.array([self.lPeptide().numberOfMatchingIons(x) + self.rPeptide().numberOfMatchingIons(x) for x in self._assignedSpectra])
  def setCheckSum(self, chksum):
    """sets the checksum for this xlink
    it should be calculated only using peptide sequences
    and their modifications
    """
    self._checksum = chksum
  def checksum(self):
    """return the checksum
    """
  
    return self._checksum
  def setPValue(self, p, sp):
    """set the p-value for match of this
    xlink to a given spectrum
    """
    
    self._pValues[sp.checksum()] = p
  def pValue(self, sp = None):
    """return the pvalue list. If a spectrum is given
    return only the corresponding p-value
    """
    if not sp:
      return self._pValues.values()
    return self._pValues[sp.checksum()]
  def setMatchScore(self, score, sp):
    """as for p-values but intensity scores
    """
    self._MatchScores[sp.checksum()] = score
  def evalMatchScores(self):
    """evaluates scores based on peptides intensity
    scores. At the same time it updates the bestmatch for every
    spectrum analyzed
    """
    for sp in self.spectra():
      score = np.mean(self.lPeptide().matchScore(sp) + self.rPeptide().matchScore(sp))
      self.setMatchScore(score , sp)
  def matchScore(self, sp = None):
    """ return the score"""  
    if not sp:
      return self._MatchScores.values()
    return self._MatchScores[sp.checksum()]
  def evalPositionCheckSum(self):
    """calculate the checksum for this xlink only based
    on mapping positions
    """
    t = set()
    lpos = self.lPeptide().positionsInProtein()
    rpos = self.rPeptide().positionsInProtein()
    lprotacc = self.lPeptide().motherProtein().accession()
    rprotacc = self.rPeptide().motherProtein().accession()
    for l in lpos:
      for r in rpos:
        m = [lprotacc + ':' + str(l + self.lXPos()), rprotacc + ':' + str(r + self.rXPos())]
        m.sort()
        c = md5(';'.join(m)).hexdigest()
        t.add(c)
    self._poschksum = list(t)
  def positionCheckSum(self, lpep = None, rpep = None, lxpos = None, rxpos = None):
    """returns the position checksums. If peptides and positions are given
    it only returns a single checksum for passed parameters
    """
    if lpep and rpep:
      lpos = lpep.positionsInProtein()
      rpos = rpep.positionsInProtein()
      lacc = lpep.motherProtein().accession()
      racc = rpep.motherProtein().accession()
      if lxpos and rxpos:
        m = [lacc + ':' + str(lxpos), racc + ':' + str(rxpos)]
        m.sort()
        c = md5(';'.join(m)).hexdigest()
      else:
        t = set()
        for l in lpos:
          for m in rpos:
            m = [lacc + ':' + str(l + self.lXPos()), racc + ':' + str(r + self.rXPos())]
            m.sort()
            c = md5(';'.join(m)).hexdigest()
            t.add(c)
        c = list(t)
      return c
    return self._poschksum
  def guessX(self, spectrum, options):
    """ Atomic function to get the best X-position for 
    both peptides given a certain spectrum
    """
    connections = XLinkers.get(self._xlinker, 'links')
    aaList = Aminoacids.sections()
    lpos2test = []
    rpos2test = []
    lscore = 1.0
    rscore = 1.0
    
    for px in connections.split(','):
      if px in aaList:
        #left peptide
        if self.lPeptide().isCTerminal(): 
          lseq = self.lPeptide().sequence()
        else:
          lseq = self.lPeptide().sequence()[:-1]
        for pos in re.finditer(px, lseq):
          lpos2test.append(pos.start() + 1)
        #right peptide
        if self.rPeptide().isCTerminal(): 
          rseq = self.rPeptide().sequence()
        else:
          rseq = self.rPeptide().sequence()[:-1]
        for pos in re.finditer(px, rseq):
          rpos2test.append(pos.start() + 1)
      elif px == 'np':
        if self.lPeptide().isNTerminal(): lpos2test.append(1)
        if self.rPeptide().isNTerminal(): rpos2test.append(1)
      elif px == 'cp':
        if self.lPeptide().isCTerminal(): lpos2test.append(self.lPeptide().length())
        if self.rPeptide().isCTerminal(): rpos2test.append(self.rPeptide().length())
  
    self.setLXPos(self.lPeptide().bestPosition(lpos2test, spectrum, options))
    self.setRXPos(self.rPeptide().bestPosition(rpos2test, spectrum, options))
  def isBestMatch(self, sp):
    """ Tells if the xlink is the best match
    for a given spectrum
    """
    c1 = sp.bestMatch().checksum()
    if self.checksum() == c1:
      return True
    return False
  def setWeight(self, w):
    self._weight = w
  def weight(self):
    return self._weight


class Peptide:
  """ class for single peptides
  """
  def __init__(self, sequence = ''):
    self._sequence = sequence
    self._positions = []
    self._modifications = {}
    self._hasXMod = False
    self._motherProtein = None
    self._xMatchedIons = {}
    self._matchScore = {}
  def copy(self):
    """ return a copy of this
    Peptide
    """
    cpep = Peptide()
    cpep._sequence = self._sequence
    cpep._positions = self._positions[:]
    cpep._modifications = self._modifications.copy()
    cpep._hasXMod = self._hasXMod
    cpep._motherProtein = self._motherProtein
    cpep._xMatchedIons = self._xMatchedIons.copy()
    cpep._matchScore = self._matchScore.copy()
    return cpep
  def setSequence(self, sequence):
    """ set peptide sequence
    """
    self._sequence = sequence.upper().strip()
  def setStartPositions(self, list):
    """ set peptide start position list on 
    a protein
    """
    self._positions = list
  def addStartPosition(self, start):
    """ add a single start position to defined list
    """
    self._positions.append(int(start))
  def addModification(self, pos, mod):
    """ add a modification to modifications dictionary
    Note: overrides existing modifications in given position
    """
    #should check that pos and mod are valid
    if pos < 1 or pos > len(self._sequence):
      raise IndexError, "Index out of peptide sequence " + mod + " "  +  str(pos) + " " + self.sequence()
    else:
      self._modifications[pos] = mod
  def setFixedModifications(self, options):
    """ On the fly fixed modification setting
    method...
    """
    for modification in options["fixedMods"]:
      modtype = Modifications.get(modification, 'modtype')
      modresidues = Modifications.get(modification, 'residues').split()
      if modresidues[0] == '*': modresidues = [chr(x) for x in range(65, 91)]
      if modtype == 'aa':
        for x in range(len(self._sequence)):
          if self._sequence[x] in modresidues:
            self.addModification(x + 1, modification)
      elif modtype == 'np' and self._sequence[0] in modresidues:
        self.addModification(1, modification)
      elif modtype == 'cp' and self._sequence[-1] in modresidues:
        self.addModification(len(self._sequence), modification)
  def modifications(self):
    """ returns modification dictionary
    """
    return self._modifications
  def numberOfModifications(self, mod):
    """ returns the number of given modification
    """
    return len([x for x in self._modifications.values() if x == mod])
  def sequence(self):
    """ return peptide sequence
    """
    return self._sequence
  def startPositions(self):
    """ return peptide start positions list
    """
    return self._positions
  def length(self):
    """ return peptide sequence length
    """
    return len(self._sequence)
  def endPositions(self):
    """ return peptide end positions list
    """
    return [x + self.length() for x in self._positions]
  def isNTerminal(self):
    """ return True if peptide is N-terminal on protein
    """
    try:
      p = self._motherProtein.sequence().find(self._sequence)
      if p == 0:
        return True
    except AttributeError:
      if 0 in self._positions:
        return True
    
    return False
  def isCTerminal(self):
    """ return True if peptide is C-terminal on protein
    """
    try:
      p = self._motherProtein.sequence().find(self._sequence)
      if p == self._motherProtein.length() - len(self._sequence):
        return True
    except AttributeError:
      return False
    
    return False
  def positionsInProtein(self):
    """return all the possible position for this peptide
    withing the mother protein
    """
    return [x.start() for x in re.finditer(self._sequence, self._motherProtein.sequence())]
  def bMass(self, ionIndex = 2):
    """ return the mass of b-ion with given ionIndex
    """
    
    if ionIndex > self.length():
      ionIndex = self.length()
    
    # check if modification in 1st aminoacid (if one) is naa or np type
    # in case set initial ion mass to mod mass
    # set it to H mass otherwise
    
    try:
      if Modifications.get(self.modifications()[1], "modtype") == "np" or Modifications.get(self.modifications()[1], "modtype") ==  "naa":
        mass = Modifications.getfloat(self.modifications()[1], "mmass")
      else:
        mass = Elements.getfloat("H", "mmass")
    except KeyError:
      mass = Elements.getfloat("H", "mmass")
    
    
    # now add all aminoacid masses plus their modification mass (if one)
    for p in range(ionIndex):
      mass += Aminoacids.getfloat(self.sequence()[p], "mmass")
      
      try:
        # modifications indeces are list indeces plus one...
        mass += Modifications.getfloat(self.modifications()[p + 1], "mmass")
      except KeyError:
        continue
    
    return mass
  def yMass(self, ionIndex = 1):
    """ return the mass of y-ion with given ionIndex
    """
    if ionIndex > self.length():
      ionIndex = self.length()
      
    mass = 0.0
#    cterm = Elements.getfloat("H", "mmass") + Elements.getfloat("C", "mmass") + 2 * Elements.getfloat("O", "mmass")
    cterm = 3 * Elements.getfloat("H", "mmass") - Elements.getfloat("e", "mmass") + Elements.getfloat("O", "mmass")
    
    # y indeces counts from right to left...
    ionIndex = self.length() - ionIndex
    
    # now add all aminoacid masses plus their modification mass (if one)
    for p in range(ionIndex, self.length()):
      mass += Aminoacids.getfloat(self.sequence()[p], "mmass")
      
      try:
        # modifications indeces are list indeces plus one...
        mass += Modifications.getfloat(self.modifications()[p + 1], "mmass")
      except KeyError:
        continue
    
    # check if modification in last aminoacid (if one) is caa or cp type
    # in case set initial ion mass to mod mass
    # set it to cterm mass otherwise
    
    try:
      if Modifications.get(self.modifications()[1], "modtype") == "cp" or Modifications.get(self.modifications()[1], "modtype") ==  "caa":
        mass += Modifications.getfloat(self.modifications()[1], "mmass")
      else:
        mass += cterm
    except KeyError:
      mass += cterm

        
    return mass
  def aMass(self, ionIndex = 2):
    """ return the mass of a-ion with given ionIndex
    """
    return self.bMass(ionIndex) - 27.9949
  def setXMod(self, b):
    self._hasXMod = bool(b)
  def hasXMod(self):
    return self._hasXMod
  def molecularMass(self):
    return self.yMass(0) + self.bMass(self.length()) - 2 * Elements.getfloat("H", "mmass")
  def gravy(self):
    """ return peptide hydropathy as described in
    Kyte J., Doolittle R.F. Reference: J. Mol. Biol. 157:105-132(1982).   
    """
    return sum([Aminoacids.getfloat(x, "gravy") for x in self.sequence()])/self.length()
  def canXLink(self, xlinker):
    """ Test if this peptide can xlink
    given a certain xlinker
    """
    for xaa in XLinkers.get(xlinker, 'links').split(','):
      if xaa in self._sequence[:-1]: return True
      if xaa == 'np' and self.isNTerminal(): return True
      if xaa == 'cp' and self.isCTerminal(): return True
    return False
  def setProtein(self, Protein):
    """ set the mother protein, i.e. the protein
    this peptide belongs to
    """
    self._motherProtein = Protein
  def motherProtein(self):
    return self._motherProtein
  def matchedIons(self, sp):
    """ return the ions matching the given spectrum
    as list of binary strings, element 0 for b-ions
    element 1 for y-ions
    """
    return self._xMatchedIons[sp.checksum()]
  def matchedBIons(self, sp):
    """ return the binary string for the matching
    b-ions
    """
    return self._xMatchedIons[sp.checksum()][0]
  def matchedYIons(self, sp):
    """ return the binary string for the matching
    y-ions
    """
    return self._xMatchedIons[sp.checksum()][1]
  def numberOfMatchingIons(self, sp):
    """ return the number of matching ions, both
    b and y ions
    """
    m = np.array([self.matchedIons(sp)], dtype=np.uint8)
    return sum(np.unpackbits(m))
  def matchScore(self, sp):
    """ return the score based on the ion intensities
    """
    return self._matchScore[sp.checksum()]
  def bestPosition(self, posList, rawSpectrum, options):
    """ Guess the best position for xlink
    in a very naive way, just count which
    xlink position raises the higher number of 
    matching ions
    """
    
    scoreList = np.ones(len(posList))
    bxMatch = np.zeros(len(posList))
    yxMatch = np.zeros(len(posList))
    posList = np.array(posList)
#    spectrum = rawSpectrum.topIons(self.length() * 4,int(rawSpectrum.peakMasses()[-1]))
    spectrum = rawSpectrum.topIons(options["topIons"],options["topWindow"])
#    if self.sequence() == 'GSTEAKITEVK' and rawSpectrum.description() == 'Elution from: 55.86 to 57.25 period: f090703_009.raw experiment: 1 cycles: 1 precIntensity: 42320146.8 FinneganScanNumber: 5210':
#      pdb.set_trace()

    for pidx in np.arange(len(posList)):
      (scoreList[pidx], bxMatch[pidx], yxMatch[pidx]) = helpers.matchIons(self, spectrum, options, splitPos = posList[pidx])
    try:
      score = scoreList.max() #reduce to the max
    except ValueError:
      return -1
   
    posList = posList[scoreList == score]

    self._xMatchedIons[rawSpectrum.checksum()] = [int(bxMatch[scoreList == score][0]), int(yxMatch[scoreList == score][0])]
    self._matchScore[rawSpectrum.checksum()] = score
       
    if len(posList) == 0:
      return -1
    else:
      return posList[0]
  def checksum(self):
    """ Return a md5 checksum based on sequence and modifications
    """
    pstring = ''
    for (k, i) in self.modifications().items():
      pstring += ':'.join([str(x) for x in [k, i]])
    pstring = self.sequence() + ';' + pstring
    return md5(pstring).hexdigest()


    
class Protein:
  """ This class defines a protein sequence
  """
  def __init__(self, definitionTuple = ('','','')):
    """ initialize the protein, possibly with a tuple
    of accession, description, sequence
    as from xllib.readEntryFromNCBI()
    """
    self._accession = definitionTuple[0]
    self._description = definitionTuple[1]
    self._sequence = definitionTuple[2]
    self._peptides = []
    self._numberOfRawPeptides = 0
  def setSequence(self, sequence):
    """ set protein sequence
    """
    self._sequence = sequence.upper().strip()
  def addSequence(self, sequence):
    """ add a sequence to previous string, useful
    when reading multiline sequencese
    """
    self._sequence += sequence.upper().strip()
  def setAccession(self, accession):
    """ set protein accession number
    """
    self._accession = accession
  def setDescription(self, description):
    """ set protein description
    """
    self._description = description
  def sequence(self):
    """ return protein sequence as string
    """
    return self._sequence
  def accession(self):
    """ return protein accession
    """
    return self._accession
  def description(self):
    """ return protein description
    """
    return self._description
  def length(self):
    """ return protein sequence length
    """
    return len(self._sequence)
  def cleave(self, options):
    """ cleaves protein with given enzyme (default Trypsin)
    and missed cleavages (default 1). Doesn't return anything
    but fills the peptide list, attribute of this class
    """
    # reset peptides...
    self._peptides = []
    missed = options["missedCleavages"]
    enzyme = options["enzyme"]
    
    if missed < 0:
      missed = 0
    
    peptides = set() #peptide set, to be unique...
    # set the first element of cleave to 0
    # so we are sure that the first peptide
    # is recognized later on
    cleavePositions = set([0])
    # get the enzyme regular expression from
    # enzymes.ini
    regexp = Enzymes.get(enzyme, "expr")
    # iterate over the sequence to get all occourrences of
    # cleavage sites
    for i in range(len(self._sequence)):
      try:
        t = re.search(regexp, self._sequence[i:]).span()
        cleavePositions.add(1 + t[0]+i)
      except AttributeError:
        continue
    cleavePositions = list(cleavePositions)
    cleavePositions.sort()
    # add the last aminoacid that has been excluded...
    cleavePositions.append(len(self._sequence))

    # get the peptides as subsequences of protein sequence
    # allowing a jump according to the number of missed cleavages
    
    for m in range(missed + 1):
      for p in range(len(cleavePositions)):
        try:
          newseq = self._sequence[cleavePositions[p]:cleavePositions[1+p+m]]
        except IndexError:
          continue
        
        if newseq not in peptides:
          peptides.add(newseq)
          if len(newseq) >= options["minLength"] and len(newseq) <= options["maxLength"]:
            newpep = Peptide()
            newpep.setSequence(newseq)
            newpep.setProtein(self)
            self._peptides.append(newpep)
    self._numberOfRawPeptides = len(self._peptides)
    peptides.clear()
      
  def peptides(self):
    """ return peptide objects
    """
    return self._peptides
  def numberOfPeptides(self):
    """ return the number of all peptide instances
    with all modifications
    """
    return len(self._peptides)
  def numberOfRawPeptides(self):
    """ return the number of peptide sequences we can have
    after cleave
    """
    return self._numberOfRawPeptides
  def items(self):
    """ return an iterator over peptides
    """
    return iter(self.peptides())
  def peptideSequences(self):
    """ return peptide sequences
    """
    return [x.sequence() for x in self.peptides()]
  def gravy(self):
    """ return peptide hydropathy as described in
    Kyte J., Doolittle R.F. Reference: J. Mol. Biol. 157:105-132(1982).   
    """
    return sum([Aminoacids.getfloat(x, "gravy") for x in self.sequence()])/self.length()
  def subsequence(self, start = 1, center = 0, span = 15):
    """ returns a peptide instance that spans a number of aminoacids
    from a starting position or centered on a specific aa
    position are given as aminoacid number
    """
    t = Peptide()
    if center:
      start = center - span / 2
    
    start -= 1    # remember that python strings are indexed starting from 0
    if start < 0:
      # if we get lower than 0... start from 0
      start = 0

    end = start + span


    if end > self.length():
      end = self.length()

    t.setSequence(self.sequence()[start:end])
    return t
  def setFixedModifications(self, options):
    """ Assign fixed mod to each peptide from this protein
    """
    [x.setFixedModifications(options) for x in self._peptides]
  def setVariableModifications(self, options):
    """ Assign variable mod to each peptide in this protein
    This method adds a new peptide instance for each combination of
    variable modifications
    """
    #for each peptide for each modification create a new combination... OMG
    bufferList = []
    for pep in self._peptides:
      tempList = [copy.deepcopy(pep)]
      positions = {}
      for modification in options["varMods"]:
        modtype = Modifications.get(modification, 'modtype')
        modresidues = Modifications.get(modification, 'residues').split()
        if modresidues[0] == '*': modresidues = [chr(x) for x in range(65, 91)]
        if modtype == 'aa':
          pattern = '[' + ''.join(modresidues) + ']'
          positions[modification] = [x.start() + 1 for x in  re.finditer(pattern, pep.sequence())]
        elif modtype == 'np' and pep.sequence()[0] in modresidues:
          positions[modification] = 1
        elif modtype == 'cp' and pep.sequence()[-1] in modresidues:
          positions[modification] = pep.length()
      for modification in positions.keys():
        for pos in positions[modification]:
          for x in range(len(tempList)):
            # create a new peptide for each present in temporary list
            # and modify it according to current modification....
            # hopefully we end up with all combinations
            t = copy.deepcopy(tempList[x])
            t.addModification(pos, modification)
            if t.numberOfModifications(modification) > options["maxVariableModifications"]: continue
            tempList.append(t)
      bufferList += tempList[1:]  # skip the first as it is the duplicate of unmodified peptide
    self._peptides += bufferList

